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We present a covariant and gauge-invariant formalism suited to the study of 
second-order effects associated with higher order tensor perturbations. The analyt¬ 
ical method we have developed enables us to characterize pure second-order tensor 
perturbations about FLRW model having different kinds of equations of state. Our 
analysis of the radiation case suggests that it may be feasible to examine the CMB 
polarization arising from higher order perturbations. 


I. INTRODUCTION 

The study of cosmology is now firmly data driven thanks to the availability of large 
amount of high quality data from numerous large-scale surveys such as those from galaxies 
red-shifts, the measurements of the CMB temperature anisotropies and polarization [1]. The 
standard model of cosmology is parametrized by six values and is based on a flat universe 
that is dominated by a cosmological constant (A) and cold dark matter (CDM), with initial 
Gaussian distribution, and inflation-seeded adiabatic fluctuations. Although this model, 
broadly speaking, successfully describes all existing CMB data [2, 3], the standard model 
raises questions among which are, what is the physics of inflation? Are the initial fluctuations 
adiabatic? The finer features of the data also raises questions on the veracity of the features 
of the underlying model. These questions demand a refinement or reexamination of the 
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theoretical and data analysis considerations. To this end, it is worth considering the role 
that of nonlinear pertnrbations might have on some of the measnrements that yield the 6 
parameters. 

Stndies of small flnctnations that are thonght to have given rise to large-scale tempera- 
tnre anisotropies and polarization of the CMB are usnally treated with hrst-order relativistic 
pertnrbation theory [4-8] . Second-order pertnrbation theory is increasingly becoming nec¬ 
essary when probing scales where linear theory becomes in accnrate. This is particnlarly 
important becanse there is no clear method within the linear theory that can help determine 
when the pertnrbations have become too large for the theory to handle [9, 10] and hence 
the need for a second order relativistic pertnrbation theory. For these reasons, second order 
pertnrbation theory has received considerable attention [9-22]. 

Second-order cosmological pertnrbations dates back to [23], who extended Lifshitz’s 
linearized theory. In that work, the anthor demonstrated that the second-order density 
contrast led to increased hrst-order density contrast over time when the pertnrbation was 
not too large. It was shown that at second order gravitational waves conld be indnced 
by deformed density pertnrbations even where hrst-order pertnrbations were non tensorial. 
Recently, the anthor of [10] examined a similar ehect in the 1-1-3 formalism. Second order 
pertnrbation theory has been nsed as primary tool for considering nonlinear dynamics. For 
example, the anthors of [15] looked at second-order pertnrbations of a hat dnst FLRW 
models with a cosmological constant. These anthors considered the evolntion of second 
order pertnrbations in hat FLRW models with A 7 ^ 0 and having a dnst eqnations of 
state. They showed that these pertnrbations tend to be constants in time, in agreement 
with the cosmic no-hair conjectnre. The anthors of [16] considered the case withont a 
cosmological constant, where gange transformations at second order were introdnced and 
nsed to stndy the gange dependency of pertnrbations. Second order ehects dnring inhation 
were stndied in [ 11 ], where the prediction of the bispectrnm of perturbation from inhation 
was examined. The full relativistic treatment of second order perturbation theory has also 
been considered in [ 8 , 18, 24, 25]. In [24], the behavior of light rays in perturbed FLRW 
models is studied and the redshift between an observer and the surface of last scattering to 
second order in the metric perturbation is explicitly calculated. Here in we hnd the problem 
associated with the size of perturbations explicitly state. In particular, it is pointed out 
that there is no guarantee that second order ehects are signihcantly smaller than those 
at hrst order, given that the large length scales associated with the problem could give 
rise to large pre-factors. The implications of this on the linear theory are signihcant and 
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need to be examined, something that we partly do in this article. Meanwhile, second order 
curvature perturbations on super-Hubble scales after inflation were studied in [22, 26]. On 
the other hand, the authors of [27, 28] have shown that second order effects may lead 
to detectable non-Gaussianity in the CMB, while those of [29, 30] have considered second 
order contributions to CMB polarization in the metric approach. The authors of [30] studied 
the B-mode CMB polarization where they found that such contributions make up part of a 
contamination in the detection of the primordial tensor modes if the tensor to scalar ratio r 
is smaller than a few xl0“®. 

The literature cited above, with a few exceptions, develop and apply metric based for¬ 
malisms in their analysis of higher-order perturbations. We are interested in a 1-1-3 covariant 
and gauge invariant approach to complement these efforts. In [10], a covariant approach to 
nonlinear perturbation theory was initiated, and the formalism used to study second order 
gravitational waves sourced by hrst order density perturbations and second order density 
perturbations sourced by gravitational waves at hrst order. However, only the dust equation 
of state was considered in that work. We extend this formalism to barotropic huids and use 
the hndings to argue for a second order effect on the CMB. This paper is organized as follows: 
In section (H) we give the background to the 1-1-3 covariant and gauge Invariant approach. 
Section (HI) discusses perturbation theory in the 1-1-3 formalism. The barotropic huid is 
considered in section (VI). Section (V) gives the analysis of the system with dust equation 
of state, complete with both analytical and numerical solutions. Section (VI) then gives the 
analysis of the case with barotropic equation of state. Here too, both the analytical and 
numerical solutions are given. Section (VH) looks at second order effects on CMB, while 
the conclusion and future work is presented in section (VHI). 


II. 1+3 FORMALISM 

All variable are dehned on a model that has a FLRW geometry of curvature 1C. This 
geometry is intrinsically linked to the 4-velocity, given by the vector tangent to the 

fundamental observer world-lines such that = dx°‘/dr and u°‘Ua = —1. Based on the 4- 
velocity, we can dehned local variables that characterize the model. To achieve this, we need 
two operators, (i) the projection tensor hab{= Qab + UaUb) which projects into the tangent 
3-spaces that are orthogonal to for the case where vorticity vanishes, (ii) We also need a 
covariant derivative Da = ();«• In general, the hrst covariant derivative of the 4-veIocity Ua 
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is given by 

DljUa ^ab 4 " (^ab ~^^hab ^a'^b^ ( 1 ) 

where 0 = D°'Ua is the expansion^ aab = D(^aUb) is the shear tensor, Uab = D\^aUb] is the 
vorticity and Aa = iia = Uh'V’^Ua is the relativistic acceleration vector ( we note that the 
spatial derivative is not projected in this case). It indicates the extent to which matter 
can be moved by none-gravitational and none-inertial forces. A perfect fluid hlled FLRW 
background is characterized by aab = ^ab = = 0. The non-zero scalars in this background 

are the expansion parameter (0), the 3-Ricci curvature JC), the energy density /i(= 

TahU°'vf’) and isotropic pressure p(= \Tabh°‘^), where Tat is the energy-momentum tensor. A 
detailed account of the 1-1-3 formalism may be found in [31, 32]. Although only the special 
case of a perfect fluid hlled model is considered in this article, the formalism can be extended 
to other types of equation of states. 

III. PERTURBATION THEORY IN THE 1-^3 FORMALISM 

The framework for perturbation theory, in the 1-1-3 formalism, was developed in [32]. The 
formalism can be considered a top-down approach in the sense that one begins with the big 
picture and then breaks it down into smaller parts. To be more precise, the methodology 
requires one to begin with propagation and constraint equations for a fully perturbed model 
then linearize about a background of choice. It is standard, and we do the same in this 
analysis, to choose the FLRW model for a background. This is because our real universe 
is, at least on large scale, well described by this model. The 1-1-3 approach to perturbation 
theory is different to the standard gauge-invariant metric perturbation theory where one 
begins with variables representing a given background and then perturbs them to desired 
order before hnding the equations of motion for perturbed model. This approach may be 
thought of as a bottom-up approach in the sense that one pieces together separate parts to 
generate more complex systems and thereby rendering the original parts constituents of the 
emergent system. Our intention is not to compare the two approaches but to provide a basic 
framework against which the work presented in this article can be understood. Imagine that 
there exist several special filters and that these Liters are able to separate gauge-invariant 
perturbations according to their sizes (or order). It should be emphasized that these Liters 
^ Let .if be a typical mean distance of some fluid behavior. From a typical volume .if^, it follows that the 
rate of change of volume is 3.if.if^. Expansion is defined as the rate of change normalized by volume. In 
particular; 3.if/.if = 0. 
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are not physical but conceptual. Of course the hrst question that needs addressing is how one 
might go about constructing them. The hrst-order hlters are constructed using conditions 
outlined by the Stewart and Walker lemma [33] . This lemma requires that a quantity vanish 
in the background, be a constant scalar held or a linear combination of products of delta 
functions for it to be gauge-invariant. One can extend this to dehne second-order-filters. In 
particular one requires the hlter to hlter-out gauge-invariant quantities that are not zeroth 
or hrst-order [34]. FIG (1) gives a schematic presentation of the approach. There is an 


Perturbed 




Background 



FIG. 1. 


inherent difficulty though in this formalism, one that is not dissimilar to that in the metric 
approach, which has to do with how to decide when a particular order of perturbation is 
‘too big’ and what to do with such. 

In the present article, we will deal with an enhancement of 1-|- 3 approach to perturbations 
to allow for the analysis of perturbations that are of the second order type. We coin the word 
Quadratization to denote the act of dropping terms of order higher than 2 in our perturbation 
scheme. This would still leave a coupling of hrst and second order perturbations and require 
a means of separating them if we wish to analyze uniquely second order effects. We will 
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describe in the next section how this may be done. 

III.l. First Order Dynamics 

First order perturbations about a FLRW background are relatively easy to characterize. 
The hrst batch are the covariant derivatives Xa = Dafi, Za = and Ca = and 

are first-order gauge-invariant (FOGI) variables corresponding respectively to the spatial 
fluctuations in the energy density, expansion rate and spatial curvature of a fluid flow. It is 
worth emphasizing that these quantities are FOGI because they are co-variantly defined and 
that they vanish exactly in the background FLRW spacetime, in line with the perturbation 
scheme outlined in the previous section. The second batch of FOGI are the shear tensor 
(Tab = D{aUh), the electric part of the Weyl tensor Eab{= CcdghC^^efh^aU‘^h^bU^ where Ccdef 
is the Weyl tensor) and the magnetic part of the Weyl tensor Hab{= \^cdefh^aU‘^h^bU^)C . 
Together, these quantities represent effects on either geometry or matter. We assume that 
the interaction between the geometry and matter is given by Einstein gravitational held 
equations (EFE): 

Eab '^Edab '^abi 

where Rab is the Ricci tensor, R is the Ricci scalar and gab the metric. Note that we have 
taken the cosmological constant A = 0. Three sets of equation emanate from the EFE and 
the integrability conditions associated with it. From the Ricci identities associated with the 
vector field u°‘ (2V[aVb]M'^ = Rab^dX^^)) obtains the Raychaudhuri equation, the vorticity 
propagation equation, the shear propagation equations, the shear divergence constraint, 
vorticity divergence identity constraint and the magnetic Weyl tensor constraint. Of these, 
only the Raychaudhuri equation (2) tells us something about the background dynamics. 

0 = -^02 - ^(/i + 3p) - 2^2 + 2a;2 + ^ ^ 2 ^ ^2) 

o ^ 

where (T^ = aabO'°‘^, = oJab^^^ and = AaA°'. Since the definition of orthogonal projec¬ 

tion onto hyper-surfaces demand that ujab = 0, it follows that = (Tab(T“^ = 0. Note that 
the background-filter renders (T^ = DaA°^ = = 0 leaving 

e = -ie2-)(f> + 3p). (3) 

This is because these terms are not zeroth-order as required of background terms. It is 
also clear that a first-order-filter will not capture cr^ or as they are products of first- 
order variables. The second-order-filter will nevertheless catch DaA°‘{^ 0). In theory, this 
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means that the different orders of perturbations could contribute to the rate of expansion in 
the fully nonlinear set up, however the contributions from these terms are often insignificant 
compared to those of zeroth-order terms and are usually neglected. This has however become 
a subject of debate in recent times [35, 36]. We assume that such back-reactions on zeroth- 
order dynamics are negligible. The twice-contracted Bianchi identities guarantee that the 
total energy-momentum is conserved. The energy conservation equation gives us the other 
background propagation equation, namely 

p = -(/i + p)0. (4) 


In a general mathematical setup, we would also require the equation for p but physical 
equation of state considerations for what we would like to study demand that either p = 0 ( 
as is the case of dust equation of state) or p =w/i ( as is the case for barotropic equation of 
state). These two options imply that p is already accounted for in the propagation equation 
for p. 

Let us now consider the consequence of first-order-filters on the general propagation 
equations for first order scalars, vectors and tensors. As mentioned above, the first order 
scalar-vector-tensor (SVT) quantities with respect to the FLRW background are aab, Eat 
and Hab- We often think of these as tensors given their transformation properties. However, 
with the exception of Hab-, the rest may be split into scalar, vector and tensor parts. The 
Hab will have vector and tensor parts, as we will describe below. We use term SVT {scalar- 
vector-tensor ) to refer to such quantities. It is easy to obtain the propagation equation for 
each of these quantities. In this section, we present them in a coupled system given by the 
matrix: 


^ab 


\ -?0 

-1 

0 


^ab 


1 

Eab 

= 

-\{p + p) 

-0 

curl 


Eab 

+ 

0 

Hab 


0 

— curl 

-0 


Hab 


0 


( 5 ) 


where curl aab = Sef{aD'^abfi (the permutation tensor Sabc = u’^Vabcd is a volume element [37]. 
The shear propagation equation is obtained from Ricci identities, while the propagation 
equations for Eab and Hab come from the Bianchi identities ViaRbc]de = 0, where the Riemann 
tensor Rabcd be split into the Ricci tensor Rab and the Weyl curvature tensor Cabcd- The 1 -|- 
3 splitting of these quantities, and EFE, together with the once-contracted Bianchi identities 
give the propagation equations for Eab and Hab [38]. We note that the full nonlinear set of 
equations have been subjected to a first-order-filter to yield matrix (5). Unlike the traditional 
matrices where the entries are either scalar or complex quantities, the square matrix in this 
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case is an operator owing to the presence of the curl terms. However, a Fonrier decomposition 
of the system wonld yield a simple sqnare matrix with time depended variables. Similarly, 
the propagation eqnations for hrst-order scalar-vector (SV) qnantities are 




1 

1 

© 

/i 





1 

2 

0 




These are derived from the time derivative of Xa and Za, the commntation relation (B3) and 
Eqs.(3) and (4). One also obtains the following FOGI constraints from the above mentioned 
Ricci and Bianchi identities. 


div((Ta, F/(j, Ha) 
curl {(Tab, ^a) 


{“^Za, 0 ) 

{Hab, 0) 


DaP = -{p+p)Ua, 


(7) 

( 8 ) 

(9) 


where div(cTa) = D^aa^- We note that constraint (9) is rednndant if either dust or barotropic 
equations of state are assumed, as we will do in this article. These sets of propagation and 
constraint equations allow the study of linear order dynamics, something that has been done 
extensively. 


III.2. High-order contributions 


Our perturbation scheme relies on reducing full nonlinear equations to a set of propagation 
and constraint equations that is interpreted as perturbations about a chosen background. 
We saw in the previous section how linearization led to the hrst order perturbations about 
the FLRW. In particular, we began with full propagation equations for each quantity and 
proceeded to drop all terms that were made of products hrst order quantities. The principle 
here is that these products are of order higher than hrst order and are too small to play 
a role at this level. The interpretation is not so simple when we begin with propagation 
equations and drop terms that are higher than second-order (i.e Quadratize). This ensures 
that we only have up to second-order terms. The subtle implication here is a back-reaction 
of second order quantities on hrst-order dynamics. This is captured in the following system: 


^ab 


r -^0 

-1 

0 


^ab 


^c(a^b) 

Eab 

= 

-\{p + p) 

-0 

curl 


Eab 

+ 

+^<^c{aEb)'^ + ‘^^cd{aA‘^Hi,'^'^ 

1 

-o 


0 

— curl —0 


1 

-O 




(10) 
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The obvious additions are the quadratic products in the last column matrix that clearly 
represent second order contribution to how FOGI variables evolve. These products are 
gauge invariant [39]. The spatial fluctuation quantities satisfy the system: 


( 11 ) 


where X = — p/(/i + p). Again the column matrix demonstrates the existence of order 
coupling, as is evident by the presence of the first order shear tensor that couples to spatial 
gradients of background scalars. This coupling, although interesting, increases the level of 
difficulty in dealing with the system. The constraint equation to the above system are given 
by: 


A’/ 


--0 

-(/i + p) -(p + p)0 




crJ’Xb 


= 

-|(/i + 3p) 

- 0 |0^ + + 3p) 



— 




-|0X 

-(/i + p)X -(/i + p)X0 




_Xa,%_ 


div(a„ Ea) = + eabccr’^'^H'^ 

diy{Ha) = i-eabca^'^E^d) 
curl (^O'ab') ^a) 0) 


( 12 ) 

(13) 

(14) 


Unlike in first-order dynamics where the different types of perturbations about FLRW (with 
p=0) decouple, the different types of perturbations are coupled in second-order dynamics. 
These coupling present several challenges. The two technical challenges encountered in 
studies of perturbations beyond the first order are (1) order mixing i.e. first - second order 
perturbations mixtures, and (2) one type of perturbation sourcing another type. These are 
in addition to problems associated with gauge issues which would arise were the variables 
not defined as gauge invariant from the beginning. We would of course want to know what 
residual gauge related issues show up at second, something that has been discussed in [39]. 
One way of skirting around these some of these problems is to develop methods for extracting 
second order gauge invariant quantities such as the ones in [40] and which we outline in the 
next section. We only present second-order pure tensor extractor as an illustration. Other 
extractors can be found in [40]. 


III.3. Higher order perturbation-type extraction 


It is a well known fact that propagation equations for scalar, vector and tensor perturba¬ 
tions couple when one considers nonlinear perturbations. One would require a method for 
extracting one type of perturbations, say tensorial type, if the objective is to analyze the 
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role that such perturbations plays in higher order dynamics. This has been a long standing 
challenge and one of the reasons why higher order perturbations have not received greater 
attention, if one were to compare with the attention given to linear perturbations. There of 
course other reasons such as the difficulty of linking higher order perturbations to cosmologi¬ 
cal observations on the one hand and the belief that linear order perturbations was sufficient 
for the link between theory and observations. These sets of propagation and constraint equa¬ 
tions allow the study of linear order dynamics. We now know that with increased precision 
in cosmological measurements and data analysis, new and hitherto unexplained features are 
becoming apparent thereby necessitating the development of high order perturbation theory 
to try to explain what linear order perturbations schemes are unable to. A comprehen¬ 
sive scheme for extracting higher order perturbations in the 1-1-3 formalism was first given 
[41] and is the scheme that we will use. The reader is referred to that article for the full 
explanation, we nevertheless present an overview here. 

The scheme is based on the assumption that all variables are defined on a background 
that has FLRW geometry of curvature /C. It can be shown that all relations developed below 
hold for objects of any perturbative order m. All VT and SVT ( rank-1 and -2 tensors ) 
given here are defined to be orthogonal to u°‘. SVT tensors are symmetric and trace-free 
(PSTF) and following [42], we use angle brackets on indices as a reminder of this fact. 
Following [41], we define a conformal spatial covariant derivative acting on scalars or spatial 
tensors as Vc = ^Dc, where ^ is the scale factor and Dc is the 1-|- 3 spatial derivative 
defined above. The operator Va has the property that it commutes with the time derivative 
operator M“Va (i.e Dc[u‘^VaVai^ = M“Va[T*cW6]) for any SVT (rank-2 tensor) Yab- 

This operator allows for the definition of the irreducible parts of the spatial derivative 
of any SVT tensors, in particular the divergence, the curl and the distortion which are 
respectfully given by: 


divn..., = 

curl Yab...c ^de{a^ ^b...c) 

disWa...?) ~ ^{cVa...b)- 


e 


(15) 

(16) 
(17) 


Any SVT tensor may then be decomposed as follows: 

Yab = Sab + Vab + Tab 


— T{a'Tb)S + V(aVfe) -|- Tab, 


(18) 


where the non-local scalar part is curl-free (curliSaf, = 0), the vector part is solenoidal 
(divV = 0 divdivV = 0 and divV^ ^ 0), while the tensor part is transverse, div7^ = 0. The 
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question of how to form local scalar, vector and tensor quantities from Yah and relate them 
to the non-local split given above was answered in [41], It was also shown that when Yah 
obeys a wave equation in the form C\Yab] = Wa6, where Wab is the source and £ contains 
time derivatives and Laplacians, and any derivative operations which preserve the rank of 
Yab - i.e., curl , disdiv or divdis , or combinations thereof. One can hnd a differential that is 
capable of extracting local tensor modes. Such an extractor is precisely what was found in 
[41] and takes the form: 

^{Yab) = [-V^ + 2/C] curl (Yab) (19) 

^ we note that curl commutes with the operator in square brackets. Applying the extractor 
to our hypothetical wave equation, it is clear that 

^{C[Yab]) = C{^[Yab]) = ^iyVab), ( 20 ) 

since ^ commutes with C. It is relatively straightforward to verify that ^ (Yab) and ^ (yVab) 
are transverse showing that the extractors yields pure tensors as required. To convert the 
extraction into Fourier space, one needs to dehne tensor harmonics 
where we have two parities of orthogonal harmonics, {k‘^ + = curl QiP « + 

-iKYI^QZ’ = curl A Fourier composed form of (20) takes the form 

(A:2 + 3/C)^^^(fc2 + 2/C)F(^) = (21) 

for a unique wave number. We will elaborate this equality when we discuss the specihc case 
later in the article. Let’s now lay down the foundation for the specihc case that we are 
interested in. 


IV. KEY SECOND ORDER VARIABLES AND EQUATIONS 


Order mixing and the mixing of perturbations types mentioned in section (HI.2) present 
signihcant challenges to the analysis of second order perturbations. The utilities discussed 
above can help resolve some of these challenges but one has to be very clear about the starting 
point if any progress is to be made. In this regard, we begin with a narrowed down case, but 
one that will help illustrate how to deal with these challenges. It is entirely plausible and 
sensible to expect that one type of perturbations at say, hrst order, should be linked to the 
same type of perturbation at say second order. This is in deed true, but it also possible for 


l-V" 


21C] can also be expressed as | 


curl - 2/C 


[41] 
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first order scalar perturbations to be linked to second order vector or tensor perturbations 
i.e. one type of perturbations acting as a source for a different type of perturbations at a 
higher order of perturbations (see the discussion of the bottom-up approach). In the top- 
down approach this says that second order tensor type may be desegregated into portions 
that link to a completely different type of perturbations at the hrst order. 

We look to 1+3 splitting and the extractor to define our basic variable. First note that 
our extractor is made up of pre-factor (—+ 2/C) multiplied by curl (= .if curl , where =if 
is the scale-factor). If the pre-factor is labeled a, then the extractor takes the simplified 
form aJfcurl. It is clear that the aJf takes on a numerical value in Fourier space if the 
scale factor, the wave number and the curvature are specified. With this in mind, we can 
dehne a simplihed second order gauge invariant quantity. We note that the vanishing of 
curl is at the core of the definition of the pure scalar part in the 1+3 splitting (see for 
example Eq. (18) and the accompanying explanation). We apply the extractor to first order 
shear tensor ^[uafe]. This will yield a pure hrst order tensor because the extractor does not 
change the order of the quantity but rather extracts the pure tensor part of it. But now 
assume that only scalar perturbations are excited at hrst order then tT[aab] = 0 at this order, 
but 7 ^ 0 at second order. This, using Stewart-Walker lemma, is a gauge-invariant 

quantity at upto second order. We denote this new quantity by = tT[crab\- Since we do 
not have vectors and tensors at hrst order ( remember we only excited scalar perturbations), 
it is sufficient to use curl as a tensor extractor rather than the full This is because the 
pre-factor (—+ 2/C), which appears before curl in Eq. (19), is rank preserving while the 
whole extractor commutes with the time derivative and may be treated as a unit. For our 
purposes, it will be sufficient to use just the curl as an extractor as we explain above. In 
general though, one would need to apply the full extractor to obtain the tensor part. Since 
the general equation of state adds another level of complexity, we will restrict our discussion 
to two simple cases: that of ’dust’ (p=0) (this has been considered elsewhere [10]) and the 
barotropic equation state (p=w/i). To the best of our knowledge, the latter development is 
new and fills a gap thereby allowing a complete formalism for perfect fluids. The case for 
imperfect fluid has not been developed and is reserved for future work. 


V. THE CASE OF DUST EQUATION OF STATE 

As discussed above, we consider second order pure tensor perturbations. We first note 
that for ’dust’ EOS w=0 (=^ p = 0 and from Eq. (9) that A“=0). Now consider the 
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special case where only linear order perturbations are excited at first order. Following 
our discussion in the previous section, curl aab characterizes second order perturbations and 
curl (Jafe ( where again a is the pre-factor mentioned above and ^ the scale factor) 
represents pure second order tensor. Using our idea of filters we can group variables into 
background quantities (/r, 0 and p = 0), hrst order quantities (scalar parts of XayZajCTab 
and Eab,) and second-order quantities (curl cxaft). Let us focus on the shear tensor, from 

constrain (Eq. 14) it is clear that curl [aab] = Hab- This means that if we did not consider 

the restriction that only scalars are present at hrst-order, then Hab would be a hrst-order 
quantity that is equivalent to the part of the shear tensor that is devoid of scalars. But 
in our restricted case, this part vanishes at hrst order, since the curl sets the scalar part 
of the shear tensor to zero. It is obvious that curl [aab] has both vector and tensor parts, 
but we can obtain a pure tensor part by applying the tensor extractor to the wave equation 
involving this quantity. Let us consider the SOGI term, 

= curl cTafe. (22) 

From Eq.(lO), it can be shown that the time derivative of T,ab{= Hab) takes the form 

Sab = -0Sab - curl (E)ab, (23) 

for the case where second-order-filters are applied coupled with the restrictions adopted in 
this is section, namely scalar perturbations at hrst-order. Now taking the time derivative of 
Eq. (23) and using the commutation relation Eq. (B5) gives 

tab + curl curl Sab + |0Sab + (0^ - h)Sab = J^ab, (24) 

where 


-Tab = + a^^ Z^|e|Tb)"] - 3curl (aaiaEbfi), 

where D\e\ implies unprojected index. We know that the left hand side is gauge invariant by 
virtual of the dehnition of Sab- The relevant question to ask is if the source term is equally 
gauge invariant. In order to ascertain this, we have to demonstrate that Eab = 0 when 
Sab = 0 (both in the background and at hrst-order filtering). There are several ways to 
achieve this but we only point out one in this article. We use the Eab pact of Eq.(lO) to 
express the last term of Eab as follows 


curl Eab + 0curl Eab + lecd{aEb)'^d\v{a'^) - curl curl Sab + ^hSab + |ccd(aCCb)‘^cliv(E'') (25) 
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This new expression can now be substituted back in the source Tab to give 

Tab = -curl Eab - ©curl Eab + curl curl -D|e|-Eb)'^, (26) 

and on applying commutation relation Eq. (B5), one hnds 

Tab = -(curl Eab)' - |0curl Eab + curl curl E^,, - \^T.ab- (27) 

We know that both E^b and Tab vanish in the background given our discussions in the 
preceding sections. Now note, from Eq. (23), that E^b = 0 ^ curlE^b = 0 and so Tab 
vanishes at first-order as required, with the implication that Tab is also gauge-invariant. 
What remains now is to solve the wave equation for E^b before applying the pure tensor 
extractor. First note that we need to determine how Tab evolves, which is what is done in 
the next section. 


V.l. Propagating the source term Tab 


It is useful to express the terms in Tab in some compact formulation that make it easier 
to handle. In this regard, we define the following new compact objects [10]. 


i^iab = ecd{adiv(cr'')crb)'^, 
llj2ab = ecd(aCliv(cr'')cTb), 
tp?.ab = ecd{a^\M{T)ab)'^, 
i’Aab = ecd(aCliv(cT'')cTfe^, 


^lab ^cd(a^ -0|e|^6) ; 

^2ab ^cd(a^ -D|e|Cr^^ , 

^3ab ^cd(a^ 

^4ab ^cd(a^ ^\e\^b) 


( 28 ) 


It follows that Tab takes the simple but compact form. 


Tab — — |'02ab — 4:Q'lpiab — ^'IpSab — 4^3ab — (29) 

Taking the time derivative of each variable given in the array (28), making use of the 
commutation relation Eq.(B5) and subjecting the resulting propagation equation to the 
second order filter yields the following closed set of first order differential equations. 





1 

1 

0 

i’lab 


3 






-(!©' -1/^) 

-20 

0 

1 

i^2ab 

= 







-(|0^ - Iai) 

0 

-20 

1 

'^lab 


0 



-B0 

3 ^ J 


V’2a6 ) 

'^lab 


( 30 ) 
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It is clear that this set takes the form xjjab = and can be solved relatively easily. Note 

that the propagation of ^|J is strnctnrally similar to that of ^ and so the ^’s will yield as system 
of the form ^ab = '^^ab', where S = d'. It is easy to obtain the solntion for the differential 
system given in Eq. (30). Note that one only needs to specify the valnes of qnantities 
0 and /i in the backgronnd. For a flat FLRW backgronnd the scale-factor ^ 

0 = 2/3f and fi = 4/9f^) where t is the proper time. Substitnting these into the system, 
integrating, isolating the relevant quantities and plugging them into J^ab, Eq. (29) yields[10], 

_7 _17 

J'ab = Oiabt -|- Pabt 3 -|- 'Jabt 3 , (31) 

where («„&, (iab and 7 ab) are the coefficients determined by the initial conditions. We present 
in appendix section (A) an alternative approach of handling the source term. What remains 
now is the solution of the non-homogeneous wave equation for J^ab- 

V.2. Solutions for the dust sub-case 

From Eqs. (31, 24 ) and the commutation relation Eq. (B6), it follows that 

^ab — -|- |0Saf, -|- (0^ —/i)Sab = + [^ab^ ^ + 'Jabt 3. (32) 

As in the previous section, both p and 0 can be expressed in terms of the proper time t. We 
require a 1-1-3 harmonics decompositions of this equation, for us to be in a position to analyze 
it. The details of these decompositions can be found in [34, 37, 43, 44]. Since T,ab and J^ab 
are functions of time and space, they should ideally be written as J2abit,x) and J^abit,x). 
The harmonic decomposition of these variables into temporal and spatial dependent parts, 
for the two parities ( electric and magnetic), take the form; 

S,b(t,x) p-)^t)Qab{K,x) + f:^^\t)Qab{K,x)] , (33) 

K 

and 

J^abit, x) = [^W (t)Q,b(«:, x) + {t)Qab{t^, x)] . (34) 

K, 

It follows that the decomposition of (32) reduces to, 

+ |0S(^) + (|^ + 0^ - (35) 

for one parity, the other parity can similarly be written down. There is hidden parity switch 
due to the effect of curl that is subsumed in this equation. In particular 

CUf\Q[ib '^Qaby 
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for a flat model. It is clear that the double curl in Eq.(24) auto-corrects this parity switch. 
It should be noted that the two parities would generally couple. Our decoupled case takes 
the form 


S(0 + (36) 

when the background variables are expressed in terms of the proper time. The challenge now 
is to pick the matching wave number for both the homogeneous part of the wave equation 
and the corresponding source. In particular, we chose identical k. It is straight forward to 
show that the solutions to the homogeneous part of the wave equation (35) are, 

£<'■'= C'iW(-irj(|,ii), Ef = C2W(-TY(|,i(), (37) 

where Ci{k,) and C 2 {k.) are the constants of integration, while J(= BesselJ) and Y(= 
BesselY) are the Bessel functions of the hrst and the second kinds respectively. The general 
solutions can then be found using Green’s method as follows, 

El-* = (|, Of) + C2(n)f^Y (|, Of) + fst’ j + f e)”’ j E‘yj'l'‘>(if. 

(38) 


For completeness, we also present numerical solutions to the wave equation in hgure (2) for 
selected initial conditions. The results in this section recover those found in [10]. We now 
turn to w 7 ^ 0 sub-case, where we obtain new results. 


VI. BAROTROPIC PERFECT FLUID 

The general non-vanishing pressure( p ^ 0) introduces a signihcant level of complexity. 
Progress can however be made for the barotropic case {p = The starting point is 

similar to that presented in the previous section where ’dust’ {p = 0) was considered. Here 
too we begin with proviso that only scalar perturbations are excited at hrst-order and 
similarly, T,ab is second order and gauge invariant. In this case the background quantities 
are: p, 0,p = p(l -|- w), hrst-order quantities are the scalar parts of Za, (Jab, Bab, and 
the second-order quantity is cur\aab- The presence of pressure introduces an extra term in 
the source to the propagation equation for T,ab- In particular. 


^ab + + curl {E)ab — 


(39) 
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FIG. 2. These curves represent numerical solutions to the gravitational wave equation (35) given 
different scaling of coefficients of the source term. The initial conditions are S(k) = = 

1,S(k)^ = 0, ^ = 25. The red curve (with arrow head) represents the numerical solution for the 
homogeneous wave. The green curve (with dots) is the solution with a = 0,/3 = .08,7 = —0.03. 
The black curve is the solution curve for GW whose source has the coefficients; a = 0.1,/? = 0,7 = 
- 0 . 02 . 


where we have used = if 0^. Using the momentum conservation equation (9) and the 
constraint involving the divergence of the electric part of the Weyl tensor from (7), we can 
now express in the following form, 

(/X+p) /i(l+w) /i(l+w) 

where w is constant is a true constant and takes on numerical values. Taking the time 
derivative of (39), using the commutation relation Eq. (B5) and applying the second-order- 
filter gives 


tab + curl curl + [0^ - /i(l + 2w)]Ea6 = Fab, (41) 


^ This should not be confused with the Aab used in (A3) 
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where Fah is given by 


ecd(a[i^fe>"cliv(a'=) + 1(1 + d\M{EY + 

(13+6w)w p. 




( 42 ) 


It is necessary to express the source Eat in terms of one of either a or This will allow 
us to construct variables analogous to those in Eq. (30). We elect to eliminate the shear 
tensor. Since a appears in Eab as part of a product (e.g crE), we use Eq. (5) to express a in 
terms of E. It follows that aY and div(cr)'^ can be expressed as 


^ d 
o^b = 


divcr'^ = 


2{E^ + QEY) 

p-(l + w) 

2 (div(E)‘= + ©divE'^ 


(43) 


(44) 


/i(l + w) 

Although these expressions do not apply when w = — 1, they suffice for the case we are inter¬ 
ested in. For now, our substitutions lead to a source that is made up of products that involve 
only the electric part of the Weyl tensor. These are ecd{aEb)'^d\vE^ and ecdiaE^'^ D\e\Eb)'^ and 
their temporal derivatives. It is easily to demonstrate that the propagation of these two 
terms give rise to two separate, but structurally similar systems. The solutions are also 
structurally identical and for this reason, will only develop and analyze the part involving 
ecd{aEb)'^^\yE'^. The source Eab equals 


—20-|-12w 20 — lOw — 6w^ ^ 


/i(l + w) 


^cd{aEb\) 


/i(l + w) 


-0e, 


cd{aEb) 


div(E‘= 


As in the previous section, can now introduce compact objects that are products of hrst- 
order quantities. In particular, we dehne and make use of the following new second-order 
variables. 


Ylab = ecd{aEh)'^^\y{EY, ilab = ^cd{aE^'' E\^\Eb)^ (45) 

Y2ab = ecd(a^6)div(E)'=, ^2ab = (^cd{aE'''' T>|e|Eb)'^, (46) 

Ysab = ecd{aEb)dl\y{EY, ^Sab = (^cdiaE'''' T'lelEj, (47) 

Yiab = fJ‘^cd{aLb)'^(i'iy{EY, ^Aab = ^cd{aL^'' D\e\Eb)'^ (48) 

Ybab = fJ‘^cd{aLb)^(i\y{EY ^5ab = ^cd{aE^^ D\e\E^'^^ (49) 

Ybab ~ ^d,6cd{aLb) dlv(T) 1^606 ~ (^cd{aL Z)|e|Tfe^ , (^ 0 ) 


where Lab = Di^a^^b)- Whereas the first two terms appear in the source term Eab-, the 

remaining arise in the propagation of the different terms that make up the source. 
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VI. 1. Analyzing the source J^ab for w 7 ^ 0 and w 7 ^ — 1 


We will use a method similar to that employed in the analysis of the source term for the 
dust case. The propagation of the variables (45-50), using the relevant relations in Eq. (5) 
and applying a second-order-filter yields a system of the form 

t/’n T fin: 

where n = 1..6 and T is the coupling matrix given by 

-|0 200 

-( 1 fw)0^ -(3 + w)0 1 -4(1+ w) 


T = 


0 

0 

0 

0 


- 2(1 + 4w)02 _(^ + 2 w )0 


0 

0 

0 


0 

0 

0 


0 

-20 


0 

0 

■(1 + w) 
1 


0 

0 

0 

0 


-(l+iw )02 _(^+w )0 -4p(l + w) 
0 0 -(|-w )0 


while the residual components are give by the column matrix 

0 
0 
0 

w/iecd(a T>Vb)'^div(E)'= 
w/iecd(a T)Vb)'^div(E)'= 

-wQfieab + wfiecdiaLb)'^ D^di\y{ay + 2\Niiecd{a D‘^ab)^d\y{Ly 
We now apply the standard 1+3 harmonic decomposition to system of equations. Again we 
only present the case for one parity. In order to do this, we introduce a new notation for 
the harmonic decomposition of our new variables. For example, 

= 5^V'i(+«:')ec+aQfe>"(fc)divg^(«:'). (51) 

k=k' k=k' 

All the terms in the above system can be decomposed similarly. First notice that 


^ , r-i2 I +d I CtTP d 


\Nldecd{a D ab) div(£')'= = 


(1 + w) 


D^(E^ + eEbnd\y{E) 


and hence. 


2 w 


(1 + w) 


e,dia [^E^ + 0 Ffe"j div(F)^ = 
w/iecd(a D^ab)‘^d\y{Ey = 


2vjk^ 


\Njdecd{a D‘^abyd\y{Ly + 2\Nij,ecd{aLb}'^ D^div 


fcr)= = 


(1 + w).if2 

2w/c^ 

(l + w)^2 

4wA;2 

(l + w).i^2 


[V'2(k, k ) + fiiiK, K')]ecd{aQb)''^^'^d\yQ^'^">y 

[fi2{K, k') + fifin, K'yecdiaQb)^^{n)d\yQ‘'{Ky 

E)4 + 2fi{K, fi:')5]ecd(agb)'^(«;)divg‘'(K'). 

(52) 
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One could attempt to solve the system in Fourier space and then use this to solve the wave 
equation in Fourier space, but physically motivated considerations reduce the amount of 
algebra required to make progress and that is what we need. In this case, we only need 
the long wave length limit of the second order quantities represented by these compact 
objects. The motivation comes the inflationary and the curvaton scenarios in the early 
universe where the wavelength of cosmological perturbations responsible for seeding the 
present cosmic structures is need be much larger than the horizon scales. We similarly argue 
that the long wavelength limit, where these terms effectively vanish, will be responsible for 
seeding an enduring second order wave equation. The net effect is that the system reduces 
to the homogeneous case V’n = Note that setting both 6 and fi to their background 

values i.e. 0 = allows the homogeneous system to be easily analyzed. This will yield 

values for each of the terms given in (45-50). Substituting the solutions obtained for ^lJlab 
and %lj 2 ab into (45) shows that J^ab is equal to 

2(-2+3w) (-7+9w) (-7+15w) (17+3w) 

3(l+w) 3(l+w) 

(53) 

^(-2 + w)(-248 + 112w + 72w2), (54) 

^(-2 + w)(-253 + 110w + 75w2), (55) 

^(-2 + w)(-253 + 104w + 69w2), (56) 

%^(-2 + w)(-248 + 130w + 90w2), (57) 

^(-2 + w)(-248 + 118w + 78w2), (58) 

%^(-2 + w)(-258 + 96w + 66w2), (59) 

e constants of integration that are determined by the initial 
conditions. Now applying the standard 1-1-3 tensor harmonic decomposition to the second- 
order wave equation gives. 


2(-2+w) 


a, t (i+w) + a t (i+w) + 


where 


“2,6 = 

“ 3,6 = 
“4,6 = 

“5,6 = 
“ 6,6 = 


and where Cn^^,n = 1,2...6 an 


£<■> + + |eEW + 10^(1 - w)ei*i = 


( 60 ) 
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for each parity. We first used commutation (B6) to convert the double curl into the Laplacian 
operator —D^. The solutions to the homogeneous part are 


ll-3w 




-\(k) 


ll-3w 


5-l-3w 
6(l+w) ’ 

^o7 ’ 

(61) 

5-|-3ui 
6(l-l-w)' 


(62) 


where again J and Y are the Bessel functions of the hrst and the kinds respectively. Here 
too, the general solutions can be found using Green’s method i.e. 


^(fc) ^ ^ ^ j j (63) 

where 


■uj = J 


( (5+3u>) K f (ll+9ui) K 

\6{i+w)^ J^o J \eii+w)^^o J 


/ (5+3u>) K f\ T f (ll+9ui) K f\ 

^ 1^6(l+w)’ ^o7 \ 6(l+w) ’ ifoV • 


(64) 

Figure (3) gives the comparison of the numerical solutions for the second-order shear wave 
equation given dust w = 0 to that given radiation w = 1/3. The analytical and numerical 
solutions demonstrate that such a complex system in tractable and results interpreted with 
the context and restrictions used. But we have to ask, how relevant would such analysis 
be to the broader held of cosmology? To answer this one would have to look at aspects of 
cosmology where second order effect might play a role. We do this in the next section. . 


VII. SECOND EFFECTS AND CMB POLARIZATION 

It is known that any mechanism which produces temperature anisotropies such as per¬ 
turbations will leave nonzero polarization of the CMB [45-52]. The contribution of linear 
order perturbations to the polarization have been characterized [53] and studied in great 
detail. 

As we pointed out in the abstract and in the introduction one would like to determine 
how a mathematical characterization of second effects such as the second-order tensors 
perturbations given in the previous section relates to effects on the CMB polarization. Since 
we have given the second order tensor perturbations in the 1-1-3 formalism, it would make 
sense to consider CMB polarization from the point of view of the 1-1-3 formalism. In the 1-1-3 
covariant approach the local angular distribution of the radiation is analyzed in terms of 
the projected symmetric and trace free (PSTF) tensor-valued multi-poles, which eliminates 
the need for harmonic functions. The resulting equations allow the analysis of the evolution 
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Gravitational wave in radiation artd dust dominated regimes 


14 


12 - 


10 




- 1 


FIG. 3. The red (dotted) line is that of second-order gravitational wave amplitude S = 
in the radiation regime, while the black line that in the matter regime. It is clear that the wave 
decays much faster in the radiation regime. The initial conditions are Sg^^ = 0.00001, = 

0, k/a = 25, a = /3 = 7 = 0.1 


of anisotropies and polarization for a general pertnrbation. We first briefly review the 1+3 
decomposition techniqnes snited to the stndy of CMB anisotropies and polarization and then 
clearly point ont how the above second order pertnrbations might show np in the polarization 
measurements. For a complete description of the formalism the reader is referred to [54]. 
For a comparative metric approach the reader is referred to [30]. We follow (Challinor’s 
approach). Our description is based on the observation made by an observer co-moving the 
4-velocity We hrst note that there exist degrees of freedom in the way can be chosen 
for a general cosmological model. To eliminate these degrees of freedoms we dehne to 
coincide with the time- like eigenvector of the matter stress-energy tensor or the 4-velocity 
of some particle species. This restriction on is necessary to ensure gauge-invariance of 
the 1+3 covariant perturbation theory. 

Now consider a photon having 4-momentum, P“, ( not to be confused with the notation 
for pressure p) and energy, E ( again not to be confused with the notation for the electric 
part of the Weyl tensor). Let the direction, relative to in which the photon is propagation 
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be given by e“. It is clear that 


= P(u“ + e“). 


Let (ei)“ and ( 62 )“ be two orthonormal ( a standard consideration in the 1+3 formalism) 
pair of polarization vectors that are orthogonal to both and e“. One can dehne a pro¬ 
jection tensor T-Lab = hab + eaCb such that (ei)a = 'H\{ei)b. We set (m“, e“, ei“, 62 “) to be 
right-handed orthonormal tetrad at the observation point. Using the polarization basis vec¬ 
tor, the observer can decompose any vector field into Stokes parameters {I,Q,U,V), all 
functions of photon energy (E) and photon direction of propagation (e“). Considering the 
way Stokes parameters transform under rotation of vectors (ei)“ and ( 62 )“, one can dehne 
the polarization tensor, Pab{E,e°‘) where the only non vanishing component is given by 


Pabie^Tie.f 


1 

2 


I + Q f/ + u\ 
U-V I-Q )' 


and where i=l,2, and j=l,2. The full tensor in terms of E and e“ can be expressed in the 
form 

Pa,{E, e“) = -hPab + ^ab + ^Ve^bce'^, 

where ^is the linear polarization. I, ^ and V are functions of E and e“. The linear 
polarization tensor can be split into electric and magnetic parts i.e. 


e“) = £ab{E, e“) + Mab{E, e“) 

For linear order perturbations, scalar perturbations contribute to the electric part of 
the polarization tensor, while tensor and vector perturbations contribute to both. The 
perturbation-type mixing and order mixing we encountered when looking at second or¬ 
der perturbations will invariably have an impact on the polarization tensor. Take the 
magnetic part of the polarization tensor for example. One can express it in the form 
M 

ab + Xih where the superscripts ( 1 ) and ( 2 ) denote the order of contributing 

tensor perturbation. The exact expression and the corresponding power spectra will be 
presented in the upcoming article [55]. It is important to note that the challenge is in 
devising techniques that could segregate the two contributions given the level of noise and 
the low level of polarization in CMB. 


VIII. CONCLUSION 

The interplay between theory and observations has spurred on interest in the the held of 
cosmology. We are witnessing a period where observations is driving this whole held of study. 
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thanks to the abundance of data from various experiments. But with the abundance comes 
the need for precision, and with precision new questions. Perhaps the improving precision 
is the greatest feature of this data analysis in this held. In principle, this should allow us to 
determine the role played by nonlinear effects in cosmology, primordial or otherwise. One 
of the greatest hindrances to the study of nonlinear relativistic effects is mathematical, for 
this reason perturbation theory remains the most promising recourse. However, nonlinear 
perturbation theory is made even more complicated by gauge issues and complex equations. 
We have discussed a new second-order gauge invariant formalism suited for studying models 
with barotropic huids. Our work extends and provides alternative analytical tools to those 
presented in [10]. 

As observed in the [10], the main difficulty encountered in the 1-1-3 covariant approach 
at second-order is linked to the gauge-issue. Whereas equation (45) can be written by 
crossing off third-order terms, they are not in general integrable. This has to do with the 
fact that unlike the metric approach which solves for operators at hrst order, the covariant 
approach only solves for physical variables. In order to integrate the second-order equations 
therefore, derivative operators must operate only on variables which vanish at hrst-order 
and in the background. We have considered the the hypothetical case where only scalar 
perturbations are excited at linear. This has allowed the characterization of a second-order 
gauge invariant variables T,ab,'ip{n)ab,^{n)ab- ^ab is equivalent to the tensor part of hrst- 
order (Tafe. The 'ip(n)ab and ^{n)ab, wliich make up the source for the T,ab wave equations, 
are products of hrst-order scalars. These second-order gauge-invariant (SOGI) variables 
highlight the two diherent ways in which second-order gauge invariant variables may arise. 
From the numerical solutions to the Hab tensor wave equation, we hnd that the magnitude 
of the tensor is much greater in the radiation dominated universe than the dust dominated. 
Since tensor perturbations contribute to CMB polarization, we have provided a possible way 
of characterizing such a contribution and will build on this in a future article [55]. 


Appendix A: Alternative method for handling the source Tab for the dust case 

An alternative method for determining Tab involves taking the time derivative of (29), 
expressing the resulting equation in terms of the original variable and isolating a new source; 


^%b T 20.7^ •A-abi 


(Al) 
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where 


•^ab — y0^'01a6 + Q'ipiab + 40('02ab + 'ipSab + P,4:ab) + §©('^206 + P,3ab) ■ 

(A2) 

Now taking the time derivative of (Al), using the propagation equation each 'ips and as 
prescribed by the commutation relation (B5) and the second order filter yields: 

•^ab + 10*4.0;, = Bab-, (A3) 


where 

2 10 32 

Bab = —‘^^i'fiab + -^fiab) -^0^(3'?/’3a6 + 3'?/’2ab + 2^2a6 + ‘^fsab) “ ^0^(3'^la6 + 2^iab). 

(A4) 

It is straight forward to show that the propagation of Bab yields, 

Bab + ^QBab - 2fiAab = 0. (A5) 

The solution to the coupled system of equations given by (Al, A3, A5) reproduces (31). 

Appendix B: Commutation Relations 

The following commutation relations are satished by any scalar, /, and PSTF-tensor 
quantity. Tab, that are dehned on a FLRW background. 



(Bl) 

eabcT\T%V^^ = -TabS^^'^TAVdp, 

(B2) 

(Vo/) = Vo/-^0Vo/-cTo'VJ, 

(B3) 

V[aVb]/ = 0, 

(B4) 


(curlT^i)- = curl (t)^ - tocurlr^j + + A‘f‘ + ieVTi,'* + 

(B6) 


curl curl (r)„, = 

(B6) 


-T'^'^acafTbd + (^‘'’'(^c{aTb)d 


cd ^ 
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